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ABSTRACT 


This  research  effort  is  directed  at  determining  the  fea¬ 
sibility  of  using  structural-sensitivity  measures  as  the  basis 
for  simplifying  large  mathematical  models  of  dynamic  systems 
(e.g.,  models  of  tactical  command,  control,  and  communication 
systems).  Toward  this  end,  a  FORTRAN  program  was  written  which 
can  be  used  to  simplify  models  characterized  by  sets  of  linear 
differential  equations.  Specifically,  the  program  determines 
the  optimal  simplified  model  (i.e.,  the  set  of  coefficients 
characterizing  the  set  of  linear  differential  equations)  of 
specified  dimension  corresponding  to  the  original  linear  model 
of  higher  dimension. 

The  algorithm,  on  which  the  FORTRAN  program  is  based, 
minimizes  an  objective  function  which  is  defined  in  terms  of 
the  structural  sensitivities  of  the  state  variables  to  be  pre¬ 
served  in  the  simplified  model.  The  program  inputs  are  the  set 
of  coefficients  which  define  <1)  the  linear  differential  equa¬ 
tions  representing  the  original  model,  (2)  the  dimension  of  the 
simplified  model  to  be  determined,  and  (3)  the  set  of  parameters 
defining  the  objective  function  to  be  minimized.  The  program 
outputs  are  the  set  of  coefficients  which  define  the  optimal 
simplified  model  and  the  value  of  the  objective  function  cor¬ 
responding  to  the  optimal  simplified  model. 


„i  j  .1  i. 
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INTRODUCTION 


An  approach  is  proposed  to  simplifying  complex  mathemat¬ 
ical  models  of  dynamic  systems.  The  approach  is  based  on  the 
concept  of  structural  sensitivity. 

The  mathematical  model  to  be  simplified  is  represented 
as  a  graph  in  which  each  link  of  the  graph  represents  the  de¬ 
pendency  of  one  system  variable  on  another.  The  state  varia¬ 
bles  of  the  system  are  selected  so  as  to  include  all  the  sys¬ 
tem  variables  of  interest  which  are  also  to  be  included  in  the 
simplified  model.  The  remaining  state  variables  are  selected 
so  as  to  minimize  the  sensitivity  of  the  interesting  state  var¬ 
iables  to  the  cutting  of  all  the  links  between  the  interesting 
state  variables  and  the  remaining  state  variables.  When  this 
minimization  is  realized,  the  optimal  simplified  model,  which 
includes  all  the  interesting  state  variables,  can  be  "cut"  out 
of  the  original  system  model.  This  process  involves  finding 
that  transformation  of  the  original  state  variables  which  pre¬ 
serves  the  interesting  state  variables  while  minimizing  the 
structural  sensitivity  of  these  state  variables  with  respect  to 
the  remaining  state  variables. 

Ir.  situations  where  the  best  possible  simplified  model 
that  can  be  cut  out  of  the  original  system  model  is  not  an 
acceptable  representation  of  the  original  system,  the  dimension 
of  the  simplified  model  can  be  increased  to  improve  its  accur- 
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acy .  In  this  case,  a  transformation  must  be  found  which  not 
only  preserves  the  interesting  state  variables,  but  which  also 
partitions  the  remaining  state  variables  into  two  sets.  One 
set  of  state  variables  from  this  partition  is  added  to  the  set 
of  interesting  state  variables  to  produce  an  augmented  set  of 
interesting  state  variables;  the  simplified  model  will  now  in¬ 
clude  this  agumented  set  and  thus  be  of  higher  dimension,  and 
of  improved  accuracy.  The  other  set  of  state  variables  from 
the  partition  is  selected  so  as  to  minimize  the  structural 
sensitivity  of  the  set  of  interesting  state  variables  with 
respect  to  this  remaining  set. 
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MODEL  SIMPLIFICATION:  THE  GENERAL  AGGREGATION  PROBLEM 


The  general  aggregation  problem  is  stated  here  for  the 
case  that  the  system  of  interest  is  well  modeled  by  a  set  of 
ordinary  differential  equations.  However,  the  extention  of 
these  results  to  those  important  systems  which  are  best  mod¬ 
eled  by  discrete-event  models  is  not  trivial,  particularly 
with  respect  to  the  computation  of  sensitivities  within  dy¬ 
namic  systems.  Additional  research  is  called  for  here. 

Consider  that  the  system  to  be  studied  is  well  modeled 
by  a  set  of  ordinary  differential  equations  given  in  canonical 
state-variable  form: 

dx  =  f  (x,u),  x(t  )  -  x  t  >y  t  (1) 

dt  x  O  O  o 

where 


X  = 

*1 

.  u  = 

,  f  = 

K.1 

• 

• 

X 

n 

u 

TO 

X 

f 

xn 

For  the  case  that  n  is  very  large,  one  is  seldom  interested  in 
observing  all  the  state  variables.  In  such  cases,  one  observes 
only  a  small  subset  of  the  state  variables:  i.e., 

q  =  g(x)  (2) 

where 


■y 

y 

q  = 

q„ 

,  g  * 

g„ 

L  Pj 

pj 
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and  p  <  n.  The  function  g  is  called  the  aggregation  function. 
The  problem  of  aggregation  is  that  of  trying  to  find  a  simpler 
model  to  generate  the  observed  variables  q  than  that  provided 
by  equations  (l)  and  (2).  Ideally,  one  seeks  an  aggregated 
model  characterized  by  function  fq  such  that 

d£  =  f  (q,u),  q(tQ)  =>  g(xQ),  t  »  tQ  (3) 

dt  ” 

where 

fq  = 

In  general  there  exists  no  f^  such  that  an  aggregated 
model  can  generate  the  observed  variables  of  a  disaggregated 
model  (If  the  state  x  is  observable  through  output  q,  then  q 
cannot,  in  general,  be  generated  by  a  system  of  dimension  less 
than  n).  In  special  cases  where  an  fq  can  be  found  such  that 
equation  (3)  is  valid,  the  aggregated  model  is  said  to  be  dy¬ 
namically  exact  to  the  disaggregated  model  with  respect  to  q. 
However,  dynamic  exactness  is  so  rare  in  practical  situations 
that,  practically,  the  problem  of  aggregation  is  that  of  find¬ 
ing  a  function  fq  that  can  be  used  to  generate  an  approxima¬ 
tion  q  to  q: 

dqa 

~  =  fq<Vu)’  qa(to)  =  8(V’  t  *  Co  (4) 

'qal 
qa  "  : 

_qap 


where 


-  O  — 


Often,  the  variables  of  interest  are  so  few  in  number 
compared  to  the  dimension  of  the  disaggregated  model  (e.g.,  one 
may  be  interested  in  only  the  average  of  all  the  state  varia¬ 
bles  in  a  complex  system  having,  say,  50,000  state  variables) 
that  there  is  little  hope  of  finding  any  fq  to  generate  a  rea¬ 
sonable  approximation  to  q.  In  such  cases,  it  is  necessary  to 
increase  the  dimension  of  the  aggregated  model.  Toward  this 
end  the  aggregation  function  is  redefined: 


where 


M 

8C1(X)' 

q ,  = 

Vr+l 

— 

Vr+l00' 

: 

.^cr 

Scr(x> 

,  V 

•  _ 

• 

.^VP 

where  qc  represents  the  variables  of  interest  and  qv  represents 
the  additional  variables  to  be  included  in  the  aggregated  model 
to  increase  model  dimension  for  purposes  of  improving  the  ap¬ 
proximation.  Thus,  function  gc  is  a  fixed  function  defining 
the  variables  of  interest  and  function  gv  is  a  function  to  be 
selected  in  the  most  advantageous  manner  in  designing  the  ag¬ 
gregated  model. 

In  trading  dynamic  exactness  for  model  simplicity,  by  ac¬ 
cepting  an  approximation  to  q,  a  difficult  problem  arises. 
Namely,  one  must  have  a  basis  for  comparing  alternate  approxi¬ 
mations.  Clearly,  the  effectiveness  of  an  approximation  is 
closely  tied  to  the  use  that  the  aggregated  model  is  to  be  put 
to.  Thus,  the  criteria  that  might  be  used  in  evaluating  ag¬ 
gregated  models  to  be  used  for  estimation  and  prediction  could 
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be  significantly  different  from  the  criteria  used  when  the  mo 
dels  are  to  be  used  for  determining  controls.  The  develop¬ 
ment  of  pertinent  criteria  for  evaluating  aggregated  models 
is  an  essential  part  of  the  aggregation  problem. 


A  QUANTITATIVE  APPROACH  TO  AGGREGATION:  STRUCTURAL  SENSITIVITY 
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The  approach  taken  to  the  aggregation  problem  is  based  on 
system  structure.  Specifically,  system  models  are  represented 
by  graphs  such  as  link-node  structures  (1),  system  diagrams  (2), 
or  signal-glow  graphs  (3,  4)  in  which  certain  points  on  the 
graphs  represent  system  variables  and  the  influence  of  one  vari¬ 
able  on  another  is  denoted  by  the  existence  of  a  path  from  that 
variable  to  the  other.  Fundamental  to  this  approach  is  the  pre¬ 
mise  that  a  proposed  aggregated  model  can  be  imbe dded  within  a 
larger  system  defined  by  the  disaggregated  system  (i.e.,  the 
funciton  f  )  and  the  aggregation  function  (i.e.,  the  functions 
gc  and  gv)*  Importantly,  in  order  that  the  proposed  aggregated 
model  exactly  generate  the  variables  of  interest  qc ,  it  is  nec¬ 
essary  that  additional  variables,  say  x^,  which  are  functions  of 
the  state  variables  of  the  disaggregated  model,  be  provided  as 
special  inputs  to  the  aggregated  model.  These  relations  be¬ 
tween  the  larger-system  variables,  x^,  and  the  proposed  aggre¬ 
gated  system  variables,  gc  and  qv,  represent  connections  in  the 
system  graph.  Perfect  aggregation  is  achieved  when  the  qc  gen-- 
erated  by  the  aggregated  model  is  totally  insensitive  to  the  ex¬ 
istence  of  these  connections. 

With  such  insensitivity,  all  connections  from  the  larger 
system  can  be  literally  cut  and  the  aggregated  model  can  be  re¬ 
moved  from  the  larger  system.  This  sensitivity  of  a  system's 
variables  to  the  cutting  of  connecting  links  is  called  structural 
sensitivity.  By  introducing  a  gain  parameter  in  such  connecting 
links  it  is  possible  to  relate  structural  sensitivities  to  the 
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well-defined  parameter  sensitivities  (e.g.,  a  link  gain  equal 
to  1  implies  the  connection  exists,  and  a  link  gain  equal  to  0 
implies  the  connection  is  broken).  The  following  example  ill¬ 


ustrates  the  proposed  approach  to  aggregation. 

Consider  a  continuous  autonomous  system  that  is  well  mod¬ 
eled  by 


dx 

dt 


fx(x) 


(x  is  an  n  vector).  We  would  like  to  design  an  aggregated  model 
of  this  system  and  we  demand  that  the  aggregated  model  generate 
a  specified  set  of  outputs  qc  defined  by  a  fixed  aggregation 
function : 


( q c  is  an  r  vector).  However,  although  we  wish  to  design  an 
r-th  order  aggregated  model  in  which  qc  is  the  state,  we  are 
willing  to  increase  the  dimension  of  the  aggregated  model  by 
adding  variables  qv  to  the  state  in  the  hope  that  the  inclusion 
of  important  dynamic  modes  in  the  aggregated  model  will  lead  to 
a  better  approximation  of  q.  The  variables  qv  are  selected  by 
the  designer  as  a  function  of  the  state  variables: 


%  =  sv(x) 

( q v  is  a  p  vector).  Thus,  here  we  seek  an  aggregated  model  of 
the  form: 


dq 

— -  =  f  (q  ,cr  ) 
dt  qc v  X  V 


IT  ■  VW 
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With  proper  care  in  selecting  the  aggregation  functions  gc  and 
gv  ,  so  as  to  avoid  algebraic  dependencies  among  the  elements  of 
qc  and  qv,  one  can  think  of  qc  and  qv  as  being  a  set  or  r+p 
state  variables  in  a  new  state  description.  Denoting  the  new 
state  vector  by  x,  we  have 


where  xA  is  an  n-p-r  vector  that  is  augmented  to  qc  and  qv  to 
complete  the  state  description.  xA  is  not  unique  and  can  be 
obtained  by  an  appropriate  transformation  from  the  original 
state  description. 

XA  =  8a(x) 

Thus,  the  transformation  from  x  to  x  is  given  by 


X  = 


r*.«i 

Vx) 

8a(x) 


gt(x.) 


where  is  a  transformation  function  and,  as  such,  has  an 
inverse :  i . e . 

^  =  g'1^) 


Differentiating  x  with  respect  to  time  gives 


dx 

dt 


V  (g  (x)) 
Vx(g£(x» 
V^(g^(x); 


dx 

dt 


where  V x(gc(x) )  is  an  rxn  matrix  such  that  each  row  of  ) 
is  the  gradient,  in  the  x  space,  of  the  corresponding  element  of 
gc(x).  Thus  for  example,  the  i-th  row  of  Vx(gc(x))  is 


Vx(gci(x))  - 


8fcl 

L3X1  SX2  ’ ' * ' ■  3XnJ 


Similarly,  Vx(gv(x))  is  a  pxn  matrix  and  Vx(gy(x))  is  an  (n-p-r)xn 
•  ,  dx 

matrix.  Since  =  f^(x)  ,  we  may  write 


Vx(gc(x)) 

Vx(sv(x))  fx(x) 

Vx(gA(x}) 


And  since  x  =  g  (x)  ,  we  obtain  a  set  of  transformed  state  equa¬ 

tions  for  the  original  system: 


Vx(gv(g"i(5t))) 

vx(gA(8;1(x))) 

«« 


fx<g;1(*» 


or,  equivalently 
dqc 

IT  =  fac(qc*W 


dT  =  fqv(VW 


dT  =  WVV 


where 


f  (q 

»q  »x  )  = 

v  (c  (g.1 

qc  c 

V  A 

x  c  t 

fqv(V 

•W  = 

fxA(qc' 

■w  = 

Wh 

Figure  1  shows  the  system  diagram  for  this  transformed  system. 
Clearly,  if  gv  and  g  can  be  selected  so  that  qc  is  completely 
independent  of  x^ ,  then  perfect  aggregation  is  achieved.  This 
independence  is  equivalent  to  being  able  to  cut  the  connection 
marked  with  "X"  without  affecting  q  .  Note  that  perfect  aggre¬ 
gation  is  achieved  if  f  ( <lc  >  qv » )  is  insensitive  to  x^.  How¬ 
ever,  in  terms  of  designing  an  aggregated  model,  this  condition 

may  be  much  too  strong.  For  example,  suppose  the  system's  oper- 

dx 

ation  is  such  that  — -  =  0  (i.e.,  f  .(q  ,q  ,x.  )  =  0)  and  f  is  such 

dt  xA  c  v  A  xA 

that  ^X^UC><1V)X^)  =  °  can  be  solved  for  x^  in  terms  of  qc  and 
q  v .  In  this  case  although  f  <lc  ><3V  >  x^  )  is  a  function  of  x^, 
the  additional  relation  relating  x^  to  qc  and  qv  makes  perfect 
aggregation  possible. 

In  situations  where  no  functions  gv  and  can  be  found 
that  desensitize  qc  to  cutting  the  connections  that  bring  x^ 
to  the  fqc  and  f^v  blocks,  an  approximation  procedure  is  sug¬ 
gested  which  is  based  on  introducing  a  cutting  parameter  a. 
Figure  2  shows  the  system  diagram  of  equations  (5)  with  such 
a  cutting  parameter:  a  =  1  gives  the  original  system:  a  =  0 
cuts  the  connections.  In  this  system,  one  may  use  the  sensitiv¬ 
ity  of  qy  to  the  cutting  parameter  a  as  an  indicator  of  the  effe 
tiveness  of  aggregation  and  proceed  to  look  for  the  functions 
gv  and  that  minimize  this  sensitivity.  It  should  be  noted 
that  the  sensitivity  of  qc  to  a  not  only  depends  on  the  aggre¬ 
gation  function  gv  but  also  on  the  set  of  state  variables  selec¬ 
ted  for  x^  :  i.e.,  on  the  selection  function  g^. 

The  approach  to  the  computation  of  the  structural  sensi- 


tivities  is  straightforward.  The  outputs  qc  are  computed  with 


a  =  1.  For  large  systems,  this  computation  could  strain  the 
capacity  of  the  computer  being  used  and,  perhaps,  prove  to  be 
impractical  in  certain  cases.  In  such  situations  where  the 
limits  of  the  computer  are  being  tested,  it  is  essential  that 
efficient  computational  algorithms  be  used  (4,  5).  By  setting 
a  =  0 ,  the  smaller  proposed  aggregated  system  is  separated  from 
the  larger  system  and  the  outputs  of  this  smaller  system,  qca, 
are  computed.  The  structural  sensitivities  are  simply  the  dif¬ 
ferences  : 

Ale  = 

Aa  \a  % 


Perfect  aggregation  (i.e.,  dynamic  exactness)  is  achieved 
when  the  structural  sensitivities  are  zero  over  the  time  inter¬ 
val  of  interest.  However,  given  that  dynamic  exactness  is  gen¬ 
erally  not  possible,  one  then  seeks  the  functions  gv  and  g^  that 
in  some  way  minimize  the  structural  sensitivities.  For  example, 
a  figure  of  merit  can  be  defined  in  terms  of  the  structural  sen¬ 
sitivities.  There  are  many  possibilities  for  defining  a  figure 
of  merit.  Some  examples  are: 


|Aqc 

“i-IaET 


-(1) 


f1  Aqc  2 

M,  =  /  (-— 1 -(x))  w(t)  dx 

2  J  Aa 
0 

M3  =  /  l£Mw(T)  dT 
J0 

f,4  ■/  dT 

0  c 


L 


Figure  of  merit  M^  stresses  the  importance  of  the  final  value 
of  qc;  M2  and  M3  consider  the  accuracy  of  qc  to  be  important 
over  the  entire  time  interval  (the  weighing  function  w  allows 
the  emphasis  to  vary  over  the  time  interval);  and  M^  is  in 
terms  a  normalized  sensitivity  for  situations  in  which  one  is 
concerned  with  percent  errors.  Clearly,  the  figure  of  merit 
to  be  used  in  any  instance  depends  on  the  system  being  modeled 


LINEAR  SYSTEMS 


Although  it  is  quite  unreasonable  to  expect  that  any 
.  3 

tactical  C  system  could  be  realistically  represented  by  a 
linear  time-invariant  model,  it  is  nevertheless  useful  to 
look  at  the  aggregation  problem  for  this  special  case.  Im¬ 
portantly,  many  large  subsystems  of  tactical  CJ  systems  are 
well  modeled  by  linear  differential  equations  and  some  pro¬ 
gress  toward  obtaining  useful  aggregated  models  can  be  made 
by  aggregating  individual  subsystems  separately.  Further, 
some  of  the  ideas  set  forth  here  are  rather  simply  illustrated 
using  linear  systems  as  examples.  However,  it  must  be  noted 
that  the  assumption  of  linearity  gives  rise  to  significant 
simplifications  that  do  not  exist  for  any  other  class  of  sys¬ 
tems  . 

Consider  the  case  that  the  system  of  interest  is  well 
modeled  by  tha  set  of  linear  differential  equations,  written 
in  matrix  form: 

=  Ax  (=  f  00) 
dt  x 

where  x  is  an  n  vector  and  A  is  an  nxn  matrix  of  constants 
with  the  element  in  the  i-th  row  and  the  j-th  column  repre¬ 
sented  by  a^ j .  The  variables  of  interest,  which  are  to  be 
outputs  of  the  aggregated  model,  are  linear  combinations  of 
the  original  set  of  state  variables: 

qc  =  V  (=  gc(x» 

where  q  is  an  r  vector  (r<n)  and  G  is  an  rxn  matrix.  Vari- 
nc  c 

ables  qc  are  to  be  state  variables  of  the  aggregated  model. 
Additional  state  variables  qv,  where  q  is  a  p  vector  (p<n-r). 
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may  be  allowed  in  the  aggregated  model  to  improve  accuracy: 
qv  =  Gvx  (=  Mx)) 

where  G  is  a  pxn  matrix.  If  q  and  q  are  considered  to  be 
state  variables  in  a  transformed  coordinate  system,  an  addition¬ 
al  n-p-r  state  variables  must  be  selected,  also  as  a  linear 
combinations  of  the  original  state  variables,  to  complete  the 
transformed  state  description: 

XA  =  GAx  (=  SA(x)) 


where  x^  is  an  n-p-r  vector  and  G^  is  an  (n-p-r)xn  matrix.  Thus, 
representing  the  new  state  description  with  the  n  vector  x,  we 
have 


x  = 


Gx  (=  gt(x)) 


where  G  is  an  nxn  nonsingular  transformation  matrix  appropri¬ 
ately  constructed  from  submatrices  Gc,  Gy,  and  G^.  Differen¬ 
tiating  the  transformation  equation  x=Gx  with  respect  to  t  gives 


dx 

dt 

where 


GAG  1x 


GqvA  =  GAG 

This  can  be  written  in  expanded 
dqc  _ 

dt  "  Gqqqc  +  Gqvqv  +  GqAXA 
dqv 

d T  =  Gvqqc  +  Gwqv  +  GvAXA 
dxA 

dt  =  G  Aqqc  +  GAvqv  +  GAAxA 


form  as 


(=  fqc(qc*VXA)) 
(=  fqv(qc,qv’XA)) 
(=  fxA(qc>qv»xA)) 


(6) 


where  the  matrix  coefficients  of  qc ,  qv>  x^  are  the  appropriate 
partitions  of  GCVj^.  Figure  3  shows  the  system  diagram  corre¬ 
sponding  to  equations  (6)  with  the  cutting  parameter  a  included 
The  objective,  of  course,  is  to  find  the  transformation  subma¬ 
trices  Gv  and  G^  that  minimize  the  sensitivity  of  qc  to  the  cut 


ting  parameter  a. 


A  DIGITAL  COMPUTER  PROGRAM 


A  FORTRAN  program  was  written  which  can  be  used  in  simpli¬ 
fying  linear  systems.  The  program  determines  the  optimal  simpli¬ 
fied  model  (i.e.,  the  coefficient  matrices  Gcc  and  Gcy)  of  speci¬ 
fied  dimension  r+p  from  the  original  linear  model  of  dimension  n 
(n>r+p).  The  program  minimizes  a  figure  of  merit  which  is  de¬ 
fined  in  terms  of  the  structural  sensitivities  and  the  time  in¬ 
terval  of  interest.  The  program  inputs  are  the  original  n-di- 
mensional  model  (i.e.,  the  coefficient  matrix  A)  and  the  dimen¬ 
sion  of  the  simplified  model  to  be  developed.  Powell's  algorithm 
(6)  is  used  to  execute  the  minimization. 

The  program  consists  of  two  major  parts  (see  Figure  4). 

One  part  consists  of  the  algorithms  required  to  compute  the  value 
of  the  figure  of  merit  from  the  original  model  ari'd  the  proposed 
aggregated  model;  the  other  part  consists  of  the  optimization 
algorithm.  Figure  5  details  the  algorithm  for  the  linear  case. 
Note  that  the  A  and  Gc  matrices  are  specified  as  input  and  Gy 
and  G^  are  determined  from  the  optimization  process.  In  the 
program  written,  Powell's  method  (6)  was  used  as  the  basis  for 
the  optimization  algorithm;  Appendix  1  details  the  optimizatic .. 
algorithm  used  to  determine  the  optimal  aggregated  system  (i.e., 
the  matrices  Gv  and  G^  )  and  the  corresponding  figure  of  merit. 
Appendix  2  gives  the  definitions  of  the  important  FORTRAN  vari¬ 
ables.  A  listing  of  the  resulting  FORTRAN  program  is  given  in  the 
Appendix  3.  The  results  of  a  sample  run  are  given  in  Appendix  4. 


k 
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CONCLUSIONS 

The  research  done  thus  far  on  using  structural  sensi¬ 
tivities  as  a  basis  for  a  measure  of  effectiveness  of  aggre¬ 
gation  seems  promising.  The  FORTRAN  program  written,  which 
determines  aggregated  models  for  linear  systems,  establishes 
the  feasibility  of  the  approach,  at  least  for  an  important 
class  of  systems.  Especially  important,  insofar  as  competi¬ 
tive  systems  is  concerned,  is  that  by  using  structural  sen¬ 
sitivities  in  the  design  of  an  aggregated  model,  attention 
must  be  given  to  all  excluded  dynamic  modes;  it  is  simply  not 
sufficient  that  the  variables  of  interest  appear  to  be  re- 
sonably  approximated.  In  minimizing  structural  sensitivities 
it  is  virtually  not  possible  to  accidently  overlook  important 
system  dynamics  in  the  aggregated  model. 

Use  of  this  FORTRAN  program  on  a  variety  of  examples 
makes  it  clear  that  additional  research  is  required  to  an¬ 
swer  important  questions  concerning  both  the  effectiveness 
and  the  applicability  of  the  approach.  As  a  first  step  for 
future  research,  I  suggest  a  complete  rewriting  of  the  aggre¬ 
gation  program  either  in  a  language  such  as  APL  so  as  to  enor¬ 
mously  simplify  the  program  or,  if  again  in  FORTRAN,  on  a  sys¬ 
tem  with  efficient  linear  algebra  software  packages.  Certainly, 
having  an  experienced  programmer  write  the  program  in  a  well- 
documented  modular  form  would  be  advisable.  In  addition,  the 
optimization  algorithm  should  be  carefully  studied  for  the  pur¬ 
pose  of  improving  its  covergence  properties.  Other  optimize- 
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tion  algorithm  should  also  be  explored  to  determine  whether  the 
efficiency  of  the  search  can  be  improved. 

With  an  improved  aggregation  program,  it  will  be  possible 
to  begin  to  categorize  models  in  terms  of  whether  or  not  they 
can  be  simplified  using  this  approach.  By  characterizing  the 
properties  of  models  that  cannot  be  aggregated  using  structural 
sensitivity  measures,  insights  will  be  obtained  that  should 
allow  us  to  either  develop  methods  to  extend  the  applicability 
of  the  approach  or  to  better  define  the  limits  of  the  approach. 

Fundamental  questions  have  been  raised  which  should  be 
the  subject  of  future  research  efforts: 

1.  To  test  the  feasibility  and  effectiveness  of  the  proposed 
approach  to  aggregation  based  on  structural  sensitivities 
by  applying  it  to  model  well-known  test  systems  from  the 
literature  (e.g.,  the  examples  used  in  references  7  and 

8  ).  This  should  provide  a  comparison  of  the  proposed 

approach  to  aggregation  to  some  of  the  existing  aggregation 
methods  in  terms  of  the  effort  involved,  the  usefulness  of 
the  resulting  model,  and  the  generality  of  the  method. 

2.  To  identify  a  suitable  subsystem  of  an  actual  Air  Force 

3 

tatical  C  system  to  be  the  subject  of  a  modeling  and  sim¬ 
ulation  effort  based  on  the  proposed  approach  to  aggrega¬ 
tion. 

3.  In  designing  aggregated  models  of  large-scale  systems,  a 
minimization  must  be  carried  out  with  respect  to  a  large 
parameter  space.  Such  minimizations  can  be  difficult. 
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especially  when  there  are  many  local  minima  to  contend  with. 
Attention  should  be  given  to  trying  to  reduce  the  dimension 
of  the  parameter  space  by'  using  the  least  number  of  parame¬ 
ters  possible  in  defining  the  transformation  functions  g 
and  g^  in  this  application. 

4.  By  actually  cutting  the  connecting  links,  the  variables  x^ 
being  fed  to  the  proposed  aggregated  model  are  actually  set 
to  zero.  This  is  of  no  concern  when  dynamic  exactness  can 
be  achieved.  However,  when  the  aggregated  model  can  only 
generate  an  approximation  to  the  variables  of  interest, 

one  should  consider  the  possibility  of  introducing  bits  in¬ 
puts  to  the  aggregated  model  at  the  points  where  the  links 
have  been  cut. 

5.  Frequently,  one  may  wish  to  constrain  the  form  of  the  aggre¬ 
gated  model,  even  at  the  cost  of  having  a  deteriorated  ag¬ 
gregated  model  or  one  of  higher  dimension.  For  example,  one 
may  require  that  the  aggregated  model  be  linear  and  time- 
invariant  so  as  to  permit  analytic  studies  of  the  model  in¬ 
stead  of,  or  in  addition  to,  simulation  studies.  Methods 
for  introducing  this  model  constraint  into  the  setting  of 
structural  sensitivities  should  be  studied.  This  possi¬ 
bility  was  briefly  considered  and  the  simple  ploy  of  re¬ 
placing  the  first  of  equations  (5)  by  the  following  equa¬ 
tion  seems  promising: 

dqc 

-3 r  ■  VW  +  [Vw*a>  -  VvV 
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where  function  f  characterizes  the  constrained  aggregated 
model.  However,  further  study  is  necessary  to  determine  how 
variations  in  the  constrained  model  parameters  affect  the 
dynamics  of  the  sensitivity  function. 

6.  Considerable  effort  should  be  directed  toward  developing 
rationales  for  various  forms  of  figures  of  merits  derived 
from  structural  sensitivities.  Many  possibilities  come  to 
mind,  including  integral  forms  (with  and  without  weighting 
functions)  and  those  based  on  final  values,  and  the  impli¬ 
cations  of  each  ought  to  be  examined,  particularly  with 
respect  to  the  relationship  of  the  effectiveness  of  the 
aggregated  model  to  the  magnitude  of  the  figure  of  merit. 

7.  Since  the  computations  of  sensitivities  are  so  much  simpler 
for  static  systems  than  for  dynamic  systems,  the  possibility 
of  designing  aggregated  models  by  examining  only  the  right- 
hand  side  of  the  canonical  state  equations  should  be  care¬ 
fully  investigated. 


8.  Linear  time-invariant  systems  should  be  studied  as  an  im¬ 


portant  special  case.  Certainly,  many  important  real  sys¬ 
tems  are  modeled  as  linear  time -invariant  systems.  How¬ 
ever,  the  fact  that  linear  t ime - in var ian t  systems  yield  to 
analysis  can  be  quite  helpful  in  developing  valuable  in¬ 


sights  into  the  implications  of  structural  sensitivities. 

9.  A  system  can  be  defined  such  that  the  sensitivities  Aqc/Act 
appear  as  the  system  outputs.  Using  this  system,  the  prob¬ 
lem  of  aggregation  can  be  cast  as  a  control  problem  in 
which  the  sensitivities  can  be  considered  to  be  error  sig- 
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nals  to  be  driven  to  zero.  Such  an  approach  to  the  aggre¬ 
gation  problem  should  be  studied.  It  seems  likely  that 
certain  aspects  of  control  theory  will  prove  useful  here. 

10.  The  proposed  approach  to  aggregation  should  be  studied 
with  respect  to  the  zero-state  response  of  systems  to  clas 
sical  test  inputs  (e.g.,  unit  steps,  sinusoids,  etc.).  It 
is  clear  that  the  aggregated  models  depend  on  the  system's 
initial  state.  Yet,  in  many  systems  it  is  unlikely  that 
certain  system  state  variables  will  ever  assume  signifi¬ 
cant  values  because  of  the  large  attenuations  between  the 
input  and  the  storage  devices  associated  with  those  state 
variables.  In  such  cases,  determining  acceptable  aggre¬ 
gated  models  might  be  simplest  by  dealing  only  with  zero- 
state  input-output  responses.  For  the  linear  case,  this 
is  equivalent  to  looking  for  aggregated  transfer  functions 

11.  A  study  should  be  made  on  the  controllability  of  systems 
using  controls  derived  from  aggregated  models.  It  seems 
that  such  a  study  makes  sense  only  if  a  weaker  definition 
of  controllability  is  used  so  as  to  take  into  account  the 
extraordinary  controls  generally  necessary  before  the  neg¬ 
lected  dynamics  modes  can  significantly  affect  the  outputs 
Particular  attention  should  be  given  to  the  role  of  the 
extra  state  variables  qy  which  are  included  in  the  aggre¬ 
gated  model  only  for  accuracy.  It  may  be  desirable  to  in¬ 


clude  some  additional  state  variables  for  purposes  of  con¬ 
trollability.  ror  competitive  systems,  it  is  especially 
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important  to  determine  the  effect  that  an  opponent  can  have 
on  controllability  and  to  try  to  design  an  aggregated  model 
so  as  to  minimize  this  effect. 

12.  A  study  should  be  made  to  extend  the  results  to  discrete- 
event  systems. 

13.  A  study  should  be  made  to  examine  the  possibility  of  using 
structural  sensitivity  measures  for  decoupling  subsystems. 
Such  an  application  seems  straightforward  in  that  the  coup¬ 
ling  of  subsystems  can  be  represented  graphically  as  links 
in  the  system  graph.  Then,  all  that  is  required  is  to  find 
the  transformations  to  minimize  the  sensitivity  of  the  sub¬ 
system  variables  to  the  cutting  of  these  links  subject  to 
the  constraints  that  particular  system  variables  be  iden¬ 
tified  with  particular  subsystems. 
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FIGURE  3.  LINEAR  SYSTEM:  SYSTEM  DIAGRAM  FOR  EQUATIONS  6 
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Figure  5  Aggregation  Algorithm  for  Linear  Models 


APPENDIX  1 


Powell's  Multiparameter  Optimization  Method 

Powell's  method,  which  is  perhaps  the  most  powerful  of 
the  direct  search  methods,  has  the  property  that  it  converges 
to  the  exact  optimum  in  a  finite  number  of  steps  when  the 
function  being  optimized  is  a  quadratic  function.  The  algor¬ 
ithm  finds  the  maximum  value  of  the  function  f(x)  and  the  m- 
vector  x  giving  that  maximum.  It  is  based  on  a  series  of  one¬ 
dimensional  searches,  each  defined  by  an  m-dimensional  direc¬ 
tion  vector  p.  The  method  of  selecting  directions  of  the  one 
dimensional  searches  lies  at  the  essence  of  the  algorithm. 

The  algorithm  proceeds  as  follows: 

11  For  each  of  the  first  m  iterations,  the  direction  vectors 
p^1^,  p^)}  ...,  p^m^  are  defined  so  as  to  step  in  each  of 

the  n  coordinate  directions . 

2)  a  is  the  size  of  the  step  taken  in  direction  p.  The  size 
of  the  steps  in  the  first  n  iterations  ,...  , 

a^m^),  however,  is  not  predetermined.  One  moves  in  each 
direction  until  the  maximum  in  that  direction  is  determined 
Thus  for  each  step,  a  single  variable  search  is  conducted 
in  one  of  the  m  directions,  using  any  efficient  single¬ 
variable  search  algorithm.  In  this  program,  a  Fibonaci 


search  was  used. 


_ _ _ * 
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3)  Generate  a  new  set  of  search  direction  vectors  as  follows 


P 

P 


(2) 

(3) 


(m-1)  (m) 

P  =  P 

(m  )  (m  )  (  0  ) 

p  =  x  -  x 

4)  Define  a  new  starting  point  x^  ^  as  the  maximum  in  the  new 

( m )  . 

p  direction 

5)  Repeat  these  steps  with  the  new  x  '  and  the  new  directions 

(1)  (2)  (m) 

P  ,  P  ,  •  •  •  j  P 


Powell  has  modified  his  algorithm  to  give  better  conver¬ 
gence  when  the  initial  guess  of  x^0^  is  a  bad  one. 

1)  Define  p^1^,  p^2\...,  p^m^  in  the  m  coordinate  directions. 
Start  at  x(0) 


2)  For  k  =  l ,  2  , .  .  . m 

x(k)  =  x( k-1 )  +  a(k)p(k) 


where  a 


(k) 


is  selected  so  as  to  maximize 


_  (k-l)  „  (k ) ( 1 ) 

f  ( x  + 


a(k)(D)>  “«)[f(x(K-1)  +  a(k)p(k)).],f( 


( k  ) 
x  ) 


3)  Find  integer  j,  l^j^m,  so  that  f(x^^)  -  f(x  ^  ) 

maximum  and  define  A=f(x^^)  -  f(x^  ^ ) 


is  a 


4)  Compute 


,  (0)  ,  (m)  (0),.  Cm)  (0). 

f3  =  f(x  +2(x  -x  ))  =  f(2x  -x  ) 


f  (  x 

f  (  X 


(0) 
(  m ) 


) 

) 


Define 


5)  If  either 


f3  *  fl 


or 


(f1-2f2  +  f3)(f1-f2-A)2  >  1^A*(f1-f3)2 

(m) 

. ,  p  and  new  x 


..  ..  +  .  (1)  (2) 
use  old  directions  p  ,  p 


for  the  next  iteration. 


6)  If  neither  inequalities  in  (5)  are  true,  define 


_  .  (m)  ( 0 ) 

p  =  ( x  -x  ) , 


.  ,  _  (0)(m)  ...  (m 

look  for  new  x  =x  +ap  where  a  minimizes  f(x 


,  ,  „  (1)  (2) 

and  define  new  direction  vectors  p  ,  p 


“  4-  ~ 


APPENDIX  2 


Important  FORTRAN  Variables 

O 

A  =  an  n  vector  derived  from  the  A  matrix 

=  all  a21 • ' • anl  a12  a22 • • ,an2 • *  *  aln  a2n*,,ann 
GC  =  an  r«n  vector  derived  from  the  Gc  matrix 

Sell  § c 2 1  * " * ® crl • *  * Scln  gc2n’*‘gcrn 

GV  =  a  p*n  vector  derived  from  the  the  G  matrix 
v  v 

gvll  gv21 *  * *gvpl* ■ *svln  gc2n‘‘‘Scpn 
GD  =  a  (n-r-p)*n  vector  derived  from  the  G^  matrix 

®  All  gA2l" 'gA(n-r-p)l“  ,gAln  g  A2n  '  *  *  SA  ( n-r-p  )n 
Y  =  a  p«n  +  (n-r-p)»n  vector  derived  from  the  Gv  and  the  G^ 

matrices  (Y  is  the  vector  to  be  determined  in  the  minimi¬ 
zation  process) 

=  [  GV :  GD  ] 

2 

G  =  an  n  vector  derived  from  the  Gc,  Gv>  and  G^  matrices 

=  [gc:gv*:Gd] 

2 

GAG  =  an  n  vector  derived  from  the  GCVA  matrix 

(G  .  =  GAG_1) 
cvA 

X  =  an  n  vector  representing  the  state  of  the  original  model 
(X  =  x) 

XHAT  =  an  vector  representing  the  state  of  the  transformed  model 
( XHAT  =  2) 


■» 


» 


f 
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APPENDIX  3 


FORTRAN  Program 


COMMON  S( 10),Y(10),P( 10, 10), V( 10) ,R( 10), X( 10, 10),FUNCT( 10) 

COMMON  EO , XX , H , J  E , FO , F , F9 , C9 , F3 , D , C8 , N , J  5 , FGRAD 

COMMON /XXX/ A ( 1 00 ) , XS ( 1 0 ) , GC ( 1 0 ) , NH , R  R , PP , DELT , FMER IT , TFINAL 

COMMON /GET 1/GCC 

COMMON/TRAJEC/JOUT 

DIMENSION  GCC(IO) 

INTEGER  RR  ,PP 

JE=0 

J0UT=0 

C(150)  INPUT  CONVERGENCE  DATA (730) 

CALL  CONVDATA 

C(350)  INPUT  INITIAL  PARAMETERS (620) 

CALL  INPPARAM 

C  INPUT  SYSTEM  PARAMETERS 

CALL  SYSPAR (A,XS,GC,RR,PP, DELT , TFINAL , NN ) 

C(230)  OUTPUT  INITIAL  PARAMETERS  AND  CONVERGENCE  DATA (840) 
rflt  I  OllTPARAM 

C  ESTIMATION  OF  PARAMETERS 

^  ************************ 

Q  ************************ 

C(370)  COMPUTE  INITIAL  F(1820) 

CALL  CMPINITF 

C(390)  INITIALIZE  DIRECTION  VECTOR  PO920) 

CALL  INDIRVEC 

C(410)  STEP (2020)  IN  EACH  DIRECTION  P 

410  CALL  STEPP 

C(430)  CHECK  FOR  CONVERGENCE (2220) 

CALL  CHEKCONV 

C  C9=1  IMPLIES  CONVERGENCE 
IF(C9.EQ.1)  GO  TO  590 

C(450)  FIND  DIRECTION  GIVING  GREATEST  INCREASE  IN  F(2420) 

CALL  RAPIDF 

C(470)  COMPUTE  TEST  POINT  F(2X(N)-X(0) ) (2530) 

CALL  TESTPNT 

C( 490)  CHECK  WHETHER  DIRECTION  VECTORS  ARE  TO  BE  CHANGED(2600) 
CALL  CHKDIRVK 

C  C8=1  IMPLIES  DIRECTION  VECTORS  ARE  TO  BE  CHANGED 
IF(C8.EQ. 1 )  GO  TO  509 

C(504)  DON'T  CHANGE  DIRECTION  VECT0R5(2710) 

CALL  DIRVECOK 

C ( 507 )  REPEAT  PROCESS  WITH  NEW  STARTING  VALUES 
GO  TO  410 

C(509)  CHANGE  DIRECTION  VECTORS(2820) 

509  CALL  NEWDIREC 

C(530)  REVISE  ESTIMATES  OF  PARAMETERS( 29 1 0 ) 

CALL  REVEST 
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C(5^0)  COMPUTE  NEW  DIRECTION  VECTORS  P(2960) 

CAL.L  COMPNEWP 

C(570)  REPEAT  ESTIMATION  PROCESS  WITH  NEW  STARTING  ESTIMATES 
GO  TO  410 
590  CONTINUE 
STOP 
END 

730  SUBROUTINE  CONVDATA 

COMMON  S(10),Y(10),P(10,10),V(10),R(10),X(10,10),FUNCT(10) 

COMMON  EO,XX,H,JE,FO,F,F9,C9,F3,D,C8,N,J5,FGRAD 
C  N=NUMBER  OF  PARAMETERS 

C  EO=THE  CONVERGENCE  FACTOR  FOR  THE  MAIN  MINIMIZATION 
C  CONVERGENCE  IMPLIES  ABS(SUM(XN-XO) )  <  EO 
C  H=THE  CONVERGENCE  FACTOR  FOR  THE  ONE-DIMENSIONAL  SEARCH 
C  JE=THE  NUMBER  OF  TIMES  THE  FUNCTION  IS  COMPUTED 
READ, N,EO,H,FG RAD 
RETURN 
END 

620  SUBROUTINE  INPPARAM 

COMMON  S(10),Y(10),P(10,10),V(10),R(10),X( 10, 10) ,FUMCT( 10) 

COMMON  E0,XX,H,JE,F0,F,F9,C9,F3,D,C8,N, J5.FGRAD 
READ, (S(I), 1=1, N) 

DO  630  1=1, N 
630  X(1,I)=S(I) 

RETURN 

END 

840  SUBROUTINE  OUTPARAM 

COMMON  S( 10) , Y( 10) , PC  10, 10) ,V( 10) ,R( 10) ,X( 10. 10) ,FUNCT( 10) 

COMMON  E0.XX,H,JE,F0,F,F9,C9,F3,D,C8,N,J5,FGRAD 
WRITE (01, 842) N 

842  FORMAT ( 'THE  NUMBER  OF  PARAMETERS  IS',15) 

WRITE (01, 843 )E0 

843  FORMATC  CONVERGENCE  FACTOR  FOR  MAIN  MINIMIZATION  :',F10.7) 
WRITE(01,844)H 

844  FORMAT ( '  CONVERGENCE  FACTOR  FOR  ONE-DIMENSIONAL  SEARCH  :',F10.7) 
WRITE(01,845)FGRAD 

845  FORMAT ( '  GRADIENT  CUTOFF  FOR  ONE-DIM.  SEARCH  :’,F10.7) 
WRITE(01,850)N 

850  FORMAT ( 'THE  ',15.'  INITIAL  PARAMETER  VALUES  ARE:’) 

WRITE(01,855)(S(I),I=1,N) 

855  FORMAT(10E12.4) 

RETURN 

END 

1820  SUBROUTINE  CMPINITF 

COMMON  S(10),Y(10),P(10,10),V(10),R(10),X(10,10),FUNCT(10) 

COMMON  EO.XX.H, JE,F0,F,F9,C9,F3,D,C8,N. J5.FGRAD 
COMMOM/XXX/A( 100) ,XS( 10) ,GC(10) ,NN,RR,PP,DELT,FMERIT,TFINAL 
INTEGER  RR.PP 
DO  1850  1=1, N 
1850  Y(I)=5(I) 

CALL  FUNCTION 

FO=F 

RETURN 

END 
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4970  SUBROUTINE  FUNCTION 

COMMON  S( 10) ,Y( 10) ,P( 10,10),V(10),R(10),X(10,10) ,FUNCT(10) 
COMMON  EO,XX,H,JE,FO,F,F9,C9,F3,D,C8,N,J5,FGRAD 
COMMON /XXX/A ( 1 00 ) , XS ( 1 0 ) , GC ( 1 0 ) , NN . RR , PP , DELT , FMERIT , TFIN AL 
INTEGER  RR,PP 
JE= JE+1 

CALL  MERIT ( A , Y , GC , XS , NN , RR , P P , DELT , FMERIT , TFIN AL ) 

F=-FMERIT 

4972  CONTINUE 
RETURN 
END 

1920  SUBROUTINE  INDIRVEC 

COMMON  S( 10) ,Y( 10) ,P( 10, 10) ,V( 10) ,R( 10) ,X( 10, 10) ,FUNCT( 10) 
COMMON  EO,XX,H,JE,FO,F,F9,C9,F3,D,C8,N,J5,FGRAD 
DO  1950  1=1, N 
DO  1950  J=1,H 
1950  P(I,J)=0 

DO  2000  1=1, N 
2000  P(I,I)=1 

RETURN 
END 

2020  SUBROUTINE  STEPP 

COMMON  S(10),Y(10),P(10,10),V(10),R(10),X(10,10),FUNCT(10) 

COMMON  EO,XX,H,JE,FO,F,F9,C9,F3,D,C8,N,J5,FGRAD 
DO  2190  1=1, N 

DO  2120  J=1,N 
V(J)=P(I,J) 

IF(I.GT.1)G0  TO  2110 
R(J)=S(J) 

GO  TO  2120 
2110  R ( J ) =X( I— 1 , J) 

2120  CONTINUE 
CALL  MAX 
FUNCT(I)=F 
DO  2180  J=1,N 
2180  X(I,J)=R( J)+XX*V( J) 

2190  CONTINUE 

F9=FUNCT(N) 

RETURN 

END 

4180  SUBROUTINE  MAX 

COMMON  S(10),Y(10),P(10,10),V(10).R(10),X(10,10).FUNCT(10) 

COMMON  EO,XX,H,JE,FO.F,F9,C9,F3,D,C8,N,J5,FGRAD 

REAL  Ml, M2 

X4=0. 1 

X0=0 

4200  M1=X0-X4 
M2=X0+X4 
XX=XO 

CALL  STEPFUNC 

QO=F 

XX=M1 

CALL  STEPFUNC 
Q1=F 
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4330 


4300 

4302 


4350 


4400 

4410 


XX=M2 


CALL  STEPFUNC 
Q2=F 

IF(QO.GT.Q1)GO  TO  4400 
IF(Q0.GT.Q2)G0  TO  4480 
WRITE(01 ,4300) 

FORMAT( 'MINIMUM  DISCOVERED  IN  SUBROUTINE  MAX') 


WRITE(6,4302)Q1,Q0,Q2 
F0RMAT('Q1=’,E12.4,'  Q0=',E12.4,' 


Q2=’,E12.4) 


IF(Q2.LT.Q1 )G0  TO  4350 


X0=M2 


GO  TO  4200 
X0=M1 


GO  TO  4200 
STOP 

IF(Q0.GT.Q2)G0  TO  4550 

X4=2*X4 

M1=X0 

Q1=Q0 

X0=M2 

QO =Q2 

M2=X0+X4 

FOLD=QO 

XX=M2 

CALL  STEPFUNC 


02=F 

IF(ABSCF-FOLD) ,GT.FGRAD)GO  TO  4430 
GO  TO  4570 
4430  GO  TO  4400 
4480  X4=2*X4 
M2=X0 
U2=Q0 
X0=M1 
Q0=Q1 
M1=X0-X4 
FOLD=QO 
XX=M1 

CALL  STEPFUNC 
Q1=F 

IF(ABS(F-FOLD).GT.FGRAD)GO  TO  4485 
GO  TO  4570 
4485  GO  TO  4330 
4550  A=M1 
B=M2 

GO  TO  4580 

Q  ##******#*###»»**«*»** 

Q  #*###**  *##*#«*###*«*# 

C  ONE-DIMENSIONAL  SEARCH  AN  INTERVAL  (A.B) 
4580  ANEW=(A+B)/2-0.1«(B-A) 

BN  EW= ( A+B ) /2+0 . 1  * ( B- A ) 

4582  FORMAT (2E 15. 5) 

XX=ANEW 
CALL  STEPFUNC 
QUJEKsF 
xX=BNEW 


-37- 


CALL  STEPFUNC 
Q2NEW=F 

IF ( Q1 NEW. GT . Q2NEW)G0  TO  4585 
A=ANEW 
GO  TO  4590 
4585  B=BNEW 

4590  IF ( (B-A) ,GT.H)GO  TO  4580 
XX=(B+A)/2 
CALL  STEPFUNC 
4570  CONTINUE 
RETURN 
END 

4760  SUBROUTINE  STEPFUNC 

COMMON  S(10),Y(10),P(10,10),V(10),R(10),X(10,10),FUNCT(10) 
COMMON  E0,XX,H,JE,F0,F,F9,C9,F3,D,C8,N,J5,FGRAD 
DO  4790  1=1, N 
4790  Y(I)=R(I)+XX*V(I) 

CALL  FUNCTION 
4795  FORMAT(4E15.5) 

RETURN 

END 

2220  SUBROUTINE  CHEKCONV 

COMMON  S(10),Y(10),P(10,10),V(10),R(10),X(10,10).FUIICT(10) 

COMMON  EO,XX,H,JE,FO,F,F9,C9,F3,D,C8,N,J5,FGRAD 

COMMON/GET1/GCC 

COMMON/TRAJEC/JOUT 

DIMENSION  GCC(IO) 

IF(ABS(FO-F9).GT.FGRAD)GO  TO  2250 
WRITE (01, 2240) 

2240  FORMAT ( 'FG RAD  >  ABS(F0-F9)') 

C9=1 

GO  TO  2298 
2250  SUM=0 

DO  2260  J= 1 ,N 

2260  SUM=SUM+ABS(X(N, J)-S( J) ) 

WRITE (01, 2290) SUM, F 

2290  FORMAT(  *SUM='  ,E12.4, '  F=\E12.4) 

WRITE(01,2294)N 

2294  FORMAT ( 'THE ',15,'  PARAMETER  VALUES  ARE: ' ) 

WRITE(01 ,2295) (X(N,J),J=1,N) 

2295  FORMAT ( 10E12.4) 

IF(SUM.GT.EO)  GO  TO  2325 
C9=1 

2298  CALL  SUMMARY 
JOUT= 1 

CALL  FUNCTION 
GO  TO  2330 
2325  C9=0 
2330  RETURN 
END 

2340  SUBROUTINE  SUMMARY 

COMMON  S( 10) ,Y( 10) ,P( 10, 10) ,V( 10) ,R( 10) ,X( 10, 10) ,FUNCT( 10) 
COMMON  EO.XX.Ii , JE,F0,F.F9 ,C9 ,F3 .D.C8 .N , J5.FGRAD 
COMMON /GET 1/GCC 
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COMMON /XXX/A ( 1 00 ) , XS ( 10) ,GC(10) ,NN ,RR,PP,DELT,FMERIT,TFINAL 
DIMENSION  GCC(IO) 

INTEGER  RR.PP 
WRITE(01,2350)F 

2350  FORMAT ( 'THE  FUNCTION  VALUE  IS',E15.5) 

CALL  GETOUT(GCC,Y,N,NN,RR,PP) 

WRITEC01 ,2390) JE 

2390  FORMAT! 'THE  FUHTION  WAS  COMPUTED  ',110,'  TIMES.') 

RETURN 

END 

2420  SUBROUTINE  RAPIDF 

COMMON  S(10),Y(10),P(10,10),V(10),R(10),X(10,10),FUNCT(10) 

COMMON  E0,XX,H.JE,F0,F.F9,C9,F3,D,CS,N,J5,FGRAD 

DI=FUNCT(1)-FC 

D=DI 

J5=1 

IFCN.EQ. 1 )G0  TO  2515 
DO  2510  K=2.N 
DI =FUNCT ( K ) -FUNCT ( K- 1 ) 

IFCDI.LT. D)GO  TO  2510 

D=DI 

J5=K 

2510  CONTINUE 
2515  CONTINUE 
RETURN 
END 

2530  SUBROUTINE  TESTPNT 

COMMON  S( 10) ,Y( 10) ,P( 10, 10) ,V( 10) ,RC 10) ,X( 10, 10) ,FUNCT( 10) 
COMMON  E0,XX,H, JE,FO,F,F9,C9,F3,D,C8,N, J5,FGRAD 
DO  2560  J= 1 ,N 
2560  Y(J)=2*X(N, J)-S( J) 

CALL  FUNCTION 

F3=F 

RETURN 

END 

2600  SUBROUTINE  CHKDIRVK 

COMMON  S(10),Y(10),P(10,10),V(10),R(10),X(10,10),FUNCT(10) 
COMMON  E0,XX.H,JE.F0.F,F9,C9.F3,D,C8,N,J5,FGRAD 
F5= (F0-2*F9+F3) * (F0-F9-D) **2 
F6=(1/2)*D*(F0-F3)**2 
IFCF3.LT. FO)GO  TO  2670 
IFCF5.LT. F6)G0  TO  2690 
2670  C8=0 

GO  TO  2700 
2690  C8=1 
2700  RETURN 
END 

2710  SUBROUTINE  DIRVECOK 

COMMON  S(10).YC10),P(10.10).VC10).R(10).X(10.10).FUNCT(10) 
COMMON  EO.XX.H , JE,F0,F.F9 ,C9 ,F3,D,C8,N , J5.FGRAD 
DO  2740  J=1,N 
2740  S(J)=X(N,J) 

F0=F9 

RETURN 

END 
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2820  SUBROUTINE  NEWDIREC 

COMMON  S(10),Y(10),P(10,10),V(10).R(10),X(10,10),FUNCT(10) 
COMMON  E0,XX,H,JE,F0,F,F9,C9,F3,D,C8,N,J5,FGRAD 
DO  2860  J=1,N 
v(j)=xcn,j)-s(J) 

2SbO  K(J)=X(N,J) 

CALL  MAX 
F0=F 
RETURN 
END 

2910  SUBROUTINE  REVEST 

COMMON  S(10),Y(10),P(10,10),V(10),R(10),X(10,10) ,FUNCT( 10) 
COMMON  E0,XX,H,JE,F0,F,F9,C9,F3,D,C8,N,J5.FGRAD 
DO  2940  J=1,N 
2940  S( J) =R( J)+XX*V( J) 

RETURN 

END 

2960  SUBROUTINE  COMPNEW 

COMMON  S(10).Y(10),P(10,10),V(10),R(10),X(10,10),FUNCT(10) 
COMMON  E0,XX,H,JE,F0,F,F9,C9,F3,D,C8,N, J5.FGRAD 
IF(J5.EQ.N)GO  TO  3050 
J6=N-1 

DO  3010  I=J5, J6 

DO  3010  J= 1 ,N 
3010  P(I,J)=PCI+1,J) 

DO  3040  J=1,N 
3040  P(N, J)=V( J) 

3050  RETURN 
END 

SUBROUTINE  GETOUT(GCC , Y , N , NN , RR , PP ) 

INTEGER  RR,PP,RP 
DIMENSION  GCC( 10) ,Y( 10) 

WRITE(01,10)N 

10  FORMAT ('THE  ',15,'  '  FINAL  PARAMETER  VALUES  ARE:’) 

WRITE(01 ,12)(Y(I),I=1,N) 

12  FORMAT(  10E12.4) 

RP=RR+PP 

WRITE(01,15)RR,RR 

15  FORMAT ( 'THE  ',15,'  X  ',15,'  AGGREGATED  SYST.  COEF.  MATRIX:’) 
DO  100  1=1, RP 
J1=I 

J2=(RP~1 )*RP+I 

WRITE (01 .200) (GCC(J) , J=J1 , J2,RP) 

100  CONTINUE 

200  F0RMATO0E12.4) 

RETURN 

END 

SUBROUTINE  MERIT(A ,Y,GC,XIHIT,N ,R ,P,DELT.FMERIT.TFINAL) 
DIMENSION  A(100).X(10),Y(10).GC(10).GD(10),G(100).GA(100) 
DIMENSION  XINITC 10) ,X!1AT( 10) ,QCB( 10) .QCA( 10) 

DIMENSION  Q C( 10) ,LXX( 10) ,MXX( 10) .GV( 10) ,GAG( 1 00) ,GCC( 10) 

COMMON /GET 1/GCC 
COmON/TRAJEC/JOUT 
INTEGER  R,P 
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C  EXTRACT  GV  FROM  Y 
K=0 

IF(P.EQ.O)GO  TO  1500 
DO  1001  1=1, P 
DO  1001  J=1,N 
K=K+1 

1001  GV(K)=Y( (1-1 )*P+J) 

1500  CONTINUE 
C  EXTRACT  GD  FROM  Y 
NRP=N-R-P 

IF(NRP.EQ.O)GO  TO  2010 
K=0 

DO  2000  1=1, NRP 
K=K+1 

GD( (I-1)*N+1)=1 
NM1=N-1 

DO  2000  J=1,NM1 
K=K+1 

2000  GD(K)=Y(P*N+(I-1)*NRP+J) 

2010  CONTINUE 

C  FORM  G  FROM  GC,  GV,  AND  GD 
K=0 

DO  6000  1=1, N 
DO  3000  J=1,R 
K=K+1 

3000  G(K)=GC( (1-1 )*R+J) 

IF(P.EQ.O)GO  TO  4500 
DO  HO 00  J=1,P 
K=K+1 

4000  G(K)=GV( (T-1 )*P+J) 

4500  CONTINUE 

IF(MRP.EQ.O)GO  TO  6000 
DO  5000  J=1 ,NRP 
K=K+1 

5000  G(K)=GD( (I-1)*NRP+J) 

6000  CONTINUE 
C  COMPUTE  XHAT(O) 

CALL  GMPRD(G,XINIT,XHAT,N,N,1) 

C  EXTRACT  QCA(O)  FROM  XHAT(O) 

K=0 

DO  6555  1=1, R 
K=K+1 

6555  QCA(K)=XIiAT(K) 

T=0 

IF(JOUT.EQ.O)GO  TO  6666 
WRITE(01 ,9550)T, (QCA(I) ,QCA(I) ,1=1 ,R) 
6666  CONTINUE 
C  OBTAIN  GAG 

CALL  GMPRD(G,A,GA,N,N,N) 

CALL  MINV(G,N ,DXX,LXX,MXX) 

CALL  GMPRD(GA,G,GAG ,N ,N,N) 

C  EXTRACT  GCC  FROM  GAG 
K=0 

KRP=R+P 
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DO  7000  1=1, KRP 
DO  7000  J=1 ,KRP 
K=K+1 

7000  GCC(K)=GAG( (I-1)*N+J) 

FMERIT=0 

1=0 

8000  IF(T.GT.TFINAL)GO  TO  9999 
C  XHAT ( T+DELT ) =XHAT ( T ) +DELT* ( GAG*XH AT ( T ) ) 

CALL  GMPRD(GAG,XHAT,X ,N,N, 1 ) 

DO  9000  1=1, N 
9000  X(I)=DELT*X(I) 

CALL  GMADD ( XHAT , X , XHAT ,  N ,  1 ) 

C  EXTRACT  QC  FROM  XHAT 
K=0 

DO  9300  1=1, R 
K=K+1 

9300  QC(K)=XHAT(K) 

C  QCA ( T+DELT ) =QC A ( T ) +DELT* ( GCC  *QCA ( T ) ) 

CALL  GMPRD(GCC,QCA,QCB,R,R, 1) 

DO  9500  1=1, R 

9500  QCB ( I ) =DELT*QCB ( I ) 

CALL  GMADD (QCA, QCB, QCA, R,1) 

C  COMPUTE  FMERIT 
DO  9600  1=1, R 

FMER IT=FMERI T+DELT* ABS ( (QCA ( I ) -QC ( I ) ) /QC ( I ) ) 

9600  CONTINUE 

IF(JOUT.£Q.O)GO  TO  9555 
TD=T+DELT 

WRITE(01,9550)TD,(QC(I),QCA(I),I=1,R) 

9550  F0RMAT(F10.6,10E12.4) 

9555  CONTINUE 
T=T+DELT 
GO  TO  8000 
9999  CONTINUE 
RETURN 
END 

SUBROUTINE  SYSPAR ( A , X , GC , R , P , DELT , TFIN AL ,N ) 

DIMENSION  A( 100) ,X( 10) ,GC( 10) 

INTEGER  R,P 

READ,  N,R,P, DELT, TFIN AL 
WRITE(01,39)N 

39  FORMAT ('DIMENSION  OF  ORIGINAL  LINEAR  SYSTEM: ',15) 

NN=N*N 

WRITE(01,79)R 

79  FORMAT ( 'NUMBER  OF  STATE  VARIABLES  OF  INTEREST:  ',15) 
WRITE(01,78)P 

78  FORMAT ('NUMBER  OF  EXTRA  STATE  VARIABLES:  ',15) 

WRITE (01, 77) DELT. TFINAL 

77  FORMAT ('INTEGRATION  INTERVALS ' ,F10.5,  ’  TIME  INTERVALS' .F10.5) 
WRITE(01 ,75) 

75  FORMAT (  'ORIGINAL  SYSTEM  A-MATRIX:') 

DO  70  1=1, N 
J  1  =  1 

J2=(N-1 )*N+I 
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READ,  (A(J) ,  J=J1 ,J2.M) 

70  WRITE(01,76)(A(J),J=J1,J2,M) 

76  FORMAT  (10F12.il) 

READ, (X(I), 1=1, N) 

WRITE(01 ,55) 

55  FORMAT ( 'THE  INITIAL  STATE  IS:') 

WRITE(01 ,59) (X(I) ,1=1 ,N) 

59  FORMATC 10E12.4) 

WRITECOl ,65)R,N 

65  FORMAT ( 'THE  ',15,'  X  ’,15,’  AGGREGATION  MATRIX:') 

DO  60  1=1, R 
J1=I 

J2=(N-1)*R+I 

READ, (GC(J),J=J1,J2,R) 

60  WRITECOl, 76)(GC(J),J=J1,J2,R) 

RETURN 

END 

/INC  GMPRD ,MINV ,GMADD 
/DATA 

2, .01, .01, .000001 

0,0,0 

3,2,0, .05,1 

0,1,0 

-257.9047,-10,0 
1,1, -5 
1,1,1 

1,. 00001,. 00001 

.00001,1,. 00001 
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DIMENSION  OF  ORIGINAL  LINEAR  SYSTEM:  3 

NUMBER  OF  STATE  VARIABLES  OF  INTEREST:  2 

NUMBER  OF  EXTRA  STATE  VARIABLES:  0 

INTEGRATION  INTERVAL=  O.OIOOO  TIME  INTERVAL:  1.00000 

ORIGINAL  SYSTEM  A-MATRIX: 

0.0000  1.0000  0.0000 

-257.9048  -10.0000  0.0000 

1.0000  1.0000  -5.0000 

THE  INITIAL  STATE  IS: 

-0.1000E  01  -0.1000E  01  -0.1000E  01 
THE  2  X  3  AGGREGATION  MATRIX: 

1.0000  1.0000  1.0000 

1.0000  0.0000  0.0000 

THE  NUMBER  OF  PARAMETERS  IS  2 

CONVERGENCE  FACTOR  FOR  MAIN  MINIMIZATION  :  0.0100000 
CONVERGENCE  FACTOR  FOR  ONE-DIMENSIONAL  SEARCH  :  0.0100000 
GRADIENT  CUTOFF  FOR  ONE-DIM.  SEARCH  :  0.0000010 
THE  2  INITIAL  PARAMETER  VALUES  ARE: 

O.OOOOE  00  O.OOOOE  00 

SUM=  0.1071E  01  F=  -0.2782E  00 
THE  2  PARAMETER  VALUES  ARE: 

0.2383E-01  -0.1047E  01 
SUM:  0.1332E  00  F:  -0.5726E-01 
THE  2  PARAMETER  VALUES  ARE: 

0.1567E-01  -0.9223E  00 
SUM:  0.9421E-02  F=  -0.3173E-01 
THE  2  PARAMETER  VALUES  ARE: 

0.1692E-01  -0.9304E  00 
THE  FUNCTION  VALUE  IS  -0.31729E-01 
THE  2  FINAL  PARAMETER  VALUES  ARE: 

0.1692E-01  -0.9304E  00 

THE  2  X  2  AGGREGATED  SYST.  COEF.  MATRIX: 

-0.7946E  01  -0.2458E  03 
0.9821E  00  -0.2038E  01 

THE  FUNTION  WAS  COMPUTED  116  TIMES. 

T  Q1(T)  QA1(T)  Q2(T)  QA2 ( T ) 

(TIME  (ORIGINAL  (AGGRE-  (ORIGINAL  (AGGREGATED 

UNITS)  SYSTEM  GATED  SYSTEM  SYSTEM 

VARIABLE)  SYSTEM  VARIABLE)  VARIABLE) 

VARIABLE  ) 

0.000000  -0.3000E  01  -0.3000E  01  -0.1000E  01  -0.1000E  01 

0.010000  -0.3010E  00  -0.3037E  00  -0.1010E  01  -0.1009E  01 

0.020000  0.2208E  01  0.2201E  01  -0.9932E  00  -0.9915E  00 

0.030000  0.4476E  01  0.4463E  01  -0.9520E  00  -0.9497E  00 

0.040000  0.6463E  01  0.6442E  01  -0.8894E  00  -0.8865E  00 

0.050000  0.8137E  01  0.8109E  01  -0.8084E  00  -0.8052E  00 

0.060000  0.9479E  01  0.9444E  01  -0.7127E  00  -0.7091E  00 

0.070000  0.1048E  02  0.1044E  02  -0.6056E  00  -0.6019E  00 

0.080000  0.1114E  02  0.1109E  02  -0.4909E  00  -0.4871E  00 

0.090000  0.1146E  02  0.1140E  02  -0.3720E  00  -0.3683E  00 

0.100000  0.1146E  02  0.H40E  02  -0.2523E  00  -0.2488E  00 

0.110000  0.1118E  02  0.1 11  IE  02  -0.1350E  00  -O.1318E  00 

0.120000  0.1062E  02  0.1055E  02  -0.2299E-01  -0.2000E-01 


0.130000  0.9834E  01  0.9760E 
0.140000  0.8854E  01  0.8778E 
0.150000  0.7720E  01  0.7642E 
0.160000  0.6472E  01  0.6394E 
0.170000  0.5150E  01  0.5074E 
0.180000  0.3795E  01  0.3721E 
0.190000  0.2444E  01  0.2372E 
0.200000  0.1130E01  0.1062E 
0.210000  -0.1136E  00  -0.1788E 
0.220000  -0.1262E  01  -0.1323E 
0.230000  -0.2291E  01  -0.2348E 
0.240000  -0.3184E  01  -0.3236E 
0.250000  -0.3928E  01  -0.3976E 
0.260000  — 0 . 4 5 1 5 E  01  -0.4558E 
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